home *** CD-ROM | disk | FTP | other *** search
/ EnigmA Amiga Run 1996 March / EnigmA AMIGA RUN 05 (1996)(G.R. Edizioni)(IT)[!][issue 1996-03][Skylink CD IV].iso / earcd / gnu / gnpltsrc.lha / specfun.c < prev    next >
C/C++ Source or Header  |  1996-01-22  |  21KB  |  844 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: specfun.c,v 1.10 1995/06/16 09:31:31 drd Exp $";
  3. #endif
  4.  
  5.  
  6. /** GNUPLOT - specfun.c
  7.  *
  8.  * Copyright (C) 1986 - 1993   Thomas Williams, Colin Kelley,
  9.  *                                              Jos van der Woude
  10.  *
  11.  * Permission to use, copy, and distribute this software and its
  12.  * documentation for any purpose with or without fee is hereby granted,
  13.  * provided that the above copyright notice appear in all copies and
  14.  * that both that copyright notice and this permission notice appear
  15.  * in supporting documentation.
  16.  *
  17.  * Permission to modify the software is granted, but not the right to
  18.  * distribute the modified code.  Modifications are to be distributed
  19.  * as patches to released version.
  20.  *
  21.  * This software is provided "as is" without express or implied warranty.
  22.  *
  23.  *
  24.  * AUTHORS
  25.  *
  26.  *   Original Software:
  27.  *   Jos van der Woude, jvdwoude@hut.nl
  28.  *
  29.  * There is a mailing list for gnuplot users. Note, however, that the
  30.  * newsgroup 
  31.  *    comp.graphics.gnuplot 
  32.  * is identical to the mailing list (they
  33.  * both carry the same set of messages). We prefer that you read the
  34.  * messages through that newsgroup, to subscribing to the mailing list.
  35.  * (If you can read that newsgroup, and are already on the mailing list,
  36.  * please send a message info-gnuplot-request@dartmouth.edu, asking to be
  37.  * removed from the mailing list.)
  38.  *
  39.  * The address for mailing to list members is
  40.  *       info-gnuplot@dartmouth.edu
  41.  * and for mailing administrative requests is 
  42.  *       info-gnuplot-request@dartmouth.edu
  43.  * The mailing list for bug reports is 
  44.  *       bug-gnuplot@dartmouth.edu
  45.  * The list of those interested in beta-test versions is
  46.  *       info-gnuplot-beta@dartmouth.edu
  47.  */
  48.  
  49. #include <math.h>
  50. #include "plot.h"
  51. #include "fnproto.h"
  52.  
  53.  
  54. extern struct value stack[STACK_DEPTH];
  55. extern int s_p;
  56. extern double zero;
  57.  
  58. #define ITMAX   100
  59. #ifdef FLT_EPSILON
  60. #define MACHEPS FLT_EPSILON /* 1.0E-08 */
  61. #else
  62. #define MACHEPS 1.0E-08
  63. #endif
  64. #ifdef FLT_MIN_EXP
  65. #define MINEXP  FLT_MIN_EXP /* -88.0 */
  66. #else
  67. #define MINEXP  -88.0
  68. #endif
  69. #ifdef FLT_MAX
  70. #define OFLOW   FLT_MAX /* 1.0E+37 */
  71. #else
  72. #define OFLOW   1.0E+37
  73. #endif
  74. #ifdef FLT_MAX_10_EXP
  75. #define XBIG    FLT_MAX_10_EXP /* 2.55E+305 */
  76. #else
  77. #define XBIG    2.55E+305
  78. #endif
  79.  
  80. /*
  81.  * Mathematical constants
  82.  */
  83. #define LNPI 1.14472988584940016
  84. #define LNSQRT2PI 0.9189385332046727
  85. #ifdef PI
  86. #undef PI
  87. #endif
  88. #define PI 3.14159265358979323846
  89. #define PNT68 0.6796875
  90. #define SQRTPI 0.9189385332046727417803297
  91. #define SQRT_TWO 1.41421356237309504880168872420969809   /* JG */
  92.  
  93. #ifndef min /* GCC ST uses inline functions */
  94. #define min(a,b) ((a) < (b) ? (a) : (b))
  95. #endif
  96.  
  97. #ifndef GAMMA
  98. static int signgam = 0;
  99. #else
  100. extern int signgam; /* this is not always declared in math.h */
  101. #endif
  102.  
  103. /* Global variables, not visible outside this file */
  104. static long     Xm1 = 2147483563L;
  105. static long     Xm2 = 2147483399L;
  106. static long     Xa1 = 40014L;
  107. static long     Xa2 = 40692L;
  108.  
  109. /* Local function declarations, not visible outside this file */
  110. static double confrac __P((double a, double b, double x));
  111. static double ibeta __P((double a, double b, double x));
  112. static double igamma __P((double a, double x));
  113. static double ranf __P((double init));
  114. static double inverse_normal_func __P((double p));
  115. static double inverse_error_func __P((double p));
  116.  
  117. #ifndef GAMMA
  118. /* Provide GAMMA function for those who do not already have one */
  119. static double lngamma __P((double z));
  120. static double lgamneg __P((double z));
  121. static double lgampos __P((double z));
  122.  
  123. /**
  124.  * from statlib, Thu Jan 23 15:02:27 EST 1992 ***
  125.  *
  126.  * This file contains two algorithms for the logarithm of the gamma function.
  127.  * Algorithm AS 245 is the faster (but longer) and gives an accuracy of about
  128.  * 10-12 significant decimal digits except for small regions around X = 1 and
  129.  * X = 2, where the function goes to zero.
  130.  * The second algorithm is not part of the AS algorithms.   It is slower but
  131.  * gives 14 or more significant decimal digits accuracy, except around X = 1
  132.  * and X = 2.   The Lanczos series from which this algorithm is derived is
  133.  * interesting in that it is a convergent series approximation for the gamma
  134.  * function, whereas the familiar series due to De Moivre (and usually wrongly
  135.  * called Stirling's approximation) is only an asymptotic approximation, as
  136.  * is the true and preferable approximation due to Stirling.
  137.  * 
  138.  * Uses Lanczos-type approximation to ln(gamma) for z > 0. Reference: Lanczos,
  139.  * C. 'A precision approximation of the gamma function', J. SIAM Numer.
  140.  * Anal., B, 1, 86-96, 1964. Accuracy: About 14 significant digits except for
  141.  * small regions in the vicinity of 1 and 2.
  142.  * 
  143.  * Programmer: Alan Miller CSIRO Division of Mathematics & Statistics
  144.  * 
  145.  * Latest revision - 17 April 1988
  146.  * 
  147.  * Additions: Translated from fortran to C, code added to handle values z < 0.
  148.  * The global variable signgam contains the sign of the gamma function.
  149.  * 
  150.  * IMPORTANT: The signgam variable contains garbage until AFTER the call to
  151.  * lngamma().
  152.  * 
  153.  * Permission granted to distribute freely for non-commercial purposes only
  154.  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
  155.  */
  156.  
  157. /* Local data, not visible outside this file 
  158. static double   a[] =
  159. {
  160.     0.9999999999995183E+00,
  161.     0.6765203681218835E+03,
  162.     -.1259139216722289E+04,
  163.     0.7713234287757674E+03,
  164.     -.1766150291498386E+03,
  165.     0.1250734324009056E+02,
  166.     -.1385710331296526E+00,
  167.     0.9934937113930748E-05,
  168.     0.1659470187408462E-06,
  169. };   */
  170.  
  171. /* from Ray Toy */
  172. static double a[] = {
  173.         .99999999999980993227684700473478296744476168282198,
  174.      676.52036812188509856700919044401903816411251975244084,
  175.    -1259.13921672240287047156078755282840836424300664868028,
  176.      771.32342877765307884865282588943070775227268469602500,
  177.     -176.61502916214059906584551353999392943274507608117860,
  178.       12.50734327868690481445893685327104972970563021816420,
  179.        -.13857109526572011689554706984971501358032683492780,
  180.         .00000998436957801957085956266828104544089848531228,
  181.         .00000015056327351493115583383579667028994545044040,
  182. };
  183.  
  184. static double   lgamneg(z)
  185. double z;
  186. {
  187.     double          tmp;
  188.  
  189.     /* Use reflection formula, then call lgampos() */
  190.     tmp = sin(z * PI);
  191.  
  192.     if (fabs(tmp) < MACHEPS) {
  193.     tmp = 0.0;
  194.     } else if (tmp < 0.0) {
  195.     tmp = -tmp;
  196.         signgam = -1;
  197.     }
  198.     return LNPI - lgampos(1.0 - z) - log(tmp);
  199.  
  200. }
  201.  
  202. static double   lgampos(z)
  203. double z;
  204. {
  205.     double          sum;
  206.     double          tmp;
  207.     int             i;
  208.  
  209.     sum = a[0];
  210.     for (i = 1, tmp = z; i < 9; i++) {
  211.         sum += a[i] / tmp;
  212.     tmp++;
  213.     }
  214.  
  215.     return log(sum) + LNSQRT2PI - z - 6.5 + (z - 0.5) * log(z + 6.5);
  216. }
  217.  
  218. static double lngamma(z)
  219. double z;
  220. {
  221.     signgam = 1;
  222.  
  223.     if (z <= 0.0)
  224.     return lgamneg(z);
  225.     else
  226.     return lgampos(z);
  227. }
  228.  
  229. #define GAMMA lngamma
  230. #endif /* GAMMA */
  231.  
  232. void f_erf()
  233. {
  234. struct value a;
  235. double x;
  236.  
  237.         x = real(pop(&a));
  238. #ifndef ERF
  239.         x = x < 0.0 ? -igamma(0.5, x * x) : igamma(0.5, x * x);
  240. #else
  241.         x = erf(x);
  242. #endif
  243.         push( Gcomplex(&a,x,0.0) );
  244. }
  245.  
  246. void f_erfc()
  247. {
  248. struct value a;
  249. double x;
  250.  
  251.         x = real(pop(&a));
  252. #ifndef ERF
  253.         x = x < 0.0 ? 1.0 + igamma(0.5, x * x) : 1.0 - igamma(0.5, x * x);
  254. #else
  255.         x = erfc(x);
  256. #endif
  257.         push( Gcomplex(&a,x,0.0) );
  258. }
  259.  
  260. void f_ibeta()
  261. {
  262. struct value a;
  263. double x;
  264. double arg1;
  265. double arg2;
  266.  
  267.     x = real(pop(&a));
  268.     arg2 = real(pop(&a));
  269.     arg1 = real(pop(&a));
  270.  
  271.     x = ibeta(arg1, arg2, x);
  272.     if(x == -1.0) {
  273.         undefined = TRUE;
  274.         push( Ginteger(&a,0) );
  275.     } else
  276.         push( Gcomplex(&a,x,0.0) );
  277. }
  278.  
  279. void f_igamma()
  280. {
  281. struct value a;
  282. double x;
  283. double arg1;
  284.  
  285.     x = real(pop(&a));
  286.     arg1 = real(pop(&a));
  287.  
  288.     x = igamma(arg1,x);
  289.     if(x == -1.0) {
  290.         undefined = TRUE;
  291.         push( Ginteger(&a,0) );
  292.     } else
  293.         push( Gcomplex(&a,x,0.0) );
  294. }
  295.  
  296. void f_gamma()
  297. {
  298. register double y;
  299. struct value a;
  300.  
  301.         y = GAMMA(real(pop(&a)));
  302.     if (y > 88.0) {
  303.         undefined = TRUE;
  304.         push( Ginteger(&a,0) );
  305.     }
  306.     else
  307.         push( Gcomplex(&a,signgam * gp_exp(y),0.0) );
  308. }
  309.  
  310. void f_lgamma()
  311. {
  312. struct value a;
  313.  
  314.         push( Gcomplex(&a, GAMMA(real(pop(&a))),0.0) );
  315. }
  316.  
  317. #ifndef BADRAND
  318.  
  319. void f_rand()
  320. {
  321. struct value a;
  322.  
  323.         push( Gcomplex(&a, ranf(real(pop(&a))),0.0) );
  324. }
  325.  
  326. #else
  327.  
  328. /* Use only to observe the effect of a "bad" random number generator. */
  329. void f_rand()
  330. {
  331. struct value a;
  332.  
  333. static unsigned int y =0;
  334. unsigned int maxran = 1000;
  335.  
  336.     (void)real(pop(&a));
  337.     y = (781*y + 387) %maxran;
  338.  
  339.     push( Gcomplex(&a, (double) y /maxran,0.0) );
  340. }
  341.  
  342. #endif
  343.  
  344. /** ibeta.c
  345.  *
  346.  *   DESCRIB   Approximate the incomplete beta function Ix(a, b).
  347.  *
  348.  *                           _
  349.  *                          |(a + b)     /x  (a-1)         (b-1)
  350.  *             Ix(a, b) = -_-------_--- * |  t     * (1 - t)     dt (a,b > 0)
  351.  *                        |(a) * |(b)   /0
  352.  *
  353.  *
  354.  *
  355.  *   CALL      p = ibeta(a, b, x)
  356.  *
  357.  *             double    a    > 0
  358.  *             double    b    > 0
  359.  *             double    x    [0, 1]
  360.  *
  361.  *   WARNING   none
  362.  *
  363.  *   RETURN    double    p    [0, 1]
  364.  *                            -1.0 on error condition
  365.  *
  366.  *   XREF      lngamma()
  367.  *
  368.  *   BUGS      none
  369.  *
  370.  *   REFERENCE The continued fraction expansion as given by
  371.  *             Abramowitz and Stegun (1964) is used.
  372.  *
  373.  * Permission granted to distribute freely for non-commercial purposes only
  374.  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
  375.  */
  376.  
  377. static double          ibeta(a, b, x)
  378. double a, b, x;
  379. {
  380.     /* Test for admissibility of arguments */
  381.     if (a <= 0.0 || b <= 0.0)
  382.     return -1.0;
  383.     if (x < 0.0 || x > 1.0)
  384.     return -1.0;;
  385.  
  386.     /* If x equals 0 or 1, return x as prob */
  387.     if (x == 0.0 || x == 1.0)
  388.     return x;
  389.  
  390.     /* Swap a, b if necessarry for more efficient evaluation */
  391.     return a < x * (a + b) ? 1.0 - confrac(b, a, 1.0 - x) : confrac(a, b, x);
  392. }
  393.  
  394. static double   confrac(a, b, x)
  395. double a, b, x;
  396. {
  397.     double          Alo = 0.0;
  398.     double          Ahi;
  399.     double          Aev;
  400.     double          Aod;
  401.     double          Blo = 1.0;
  402.     double          Bhi = 1.0;
  403.     double          Bod = 1.0;
  404.     double          Bev = 1.0;
  405.     double          f;
  406.     double          fold;
  407.     double          Apb = a + b;
  408.     double          d;
  409.     int             i;
  410.     int             j;
  411.  
  412.     /* Set up continued fraction expansion evaluation. */
  413.     Ahi = gp_exp(GAMMA(Apb) + a * log(x) + b * log(1.0 - x) -
  414.               GAMMA(a + 1.0) - GAMMA(b));
  415.  
  416.     /*
  417.      * Continued fraction loop begins here. Evaluation continues until
  418.      * maximum iterations are exceeded, or convergence achieved.
  419.      */
  420.     for (i = 0, j = 1, f = Ahi; i <= ITMAX; i++, j++) {
  421.     d = a + j + i;
  422.     Aev = -(a + i) * (Apb + i) * x / d / (d - 1.0);
  423.     Aod = j * (b - j) * x / d / (d + 1.0);
  424.     Alo = Bev * Ahi + Aev * Alo;
  425.     Blo = Bev * Bhi + Aev * Blo;
  426.     Ahi = Bod * Alo + Aod * Ahi;
  427.     Bhi = Bod * Blo + Aod * Bhi;
  428.  
  429.     if (fabs(Bhi) < MACHEPS)
  430.         Bhi = 0.0;
  431.  
  432.     if (Bhi != 0.0) {
  433.         fold = f;
  434.         f = Ahi / Bhi;
  435.         if (fabs(f - fold) < fabs(f) * MACHEPS)
  436.         return f;
  437.     }
  438.     }
  439.  
  440.     return -1.0;
  441. }
  442.  
  443. /** igamma.c
  444.  *
  445.  *   DESCRIB   Approximate the incomplete gamma function P(a, x).
  446.  *
  447.  *                         1     /x  -t   (a-1)
  448.  *             P(a, x) = -_--- * |  e  * t     dt      (a > 0)
  449.  *                       |(a)   /0
  450.  *
  451.  *   CALL      p = igamma(a, x)
  452.  *
  453.  *             double    a    >  0
  454.  *             double    x    >= 0
  455.  *
  456.  *   WARNING   none
  457.  *
  458.  *   RETURN    double    p    [0, 1]
  459.  *                            -1.0 on error condition
  460.  *
  461.  *   XREF      lngamma()
  462.  *
  463.  *   BUGS      Values 0 <= x <= 1 may lead to inaccurate results.
  464.  *
  465.  *   REFERENCE ALGORITHM AS239  APPL. STATIST. (1988) VOL. 37, NO. 3
  466.  *
  467.  * Permission granted to distribute freely for non-commercial purposes only
  468.  * Copyright (c) 1992 Jos van der Woude, jvdwoude@hut.nl
  469.  */
  470.  
  471. /* Global variables, not visible outside this file */
  472. static double   pn1, pn2, pn3, pn4, pn5, pn6;
  473.  
  474. static double          igamma(a, x)
  475. double a, x;
  476. {
  477.     double          arg;
  478.     double          aa;
  479.     double          an;
  480.     double          b;
  481.     int             i;
  482.  
  483.     /* Check that we have valid values for a and x */
  484.     if (x < 0.0 || a <= 0.0)
  485.     return -1.0;
  486.  
  487.     /* Deal with special cases */
  488.     if (x == 0.0)
  489.     return 0.0;
  490.     if (x > XBIG)
  491.     return 1.0;
  492.  
  493.     /* Check value of factor arg */
  494.     arg = a * log(x) - x - GAMMA(a + 1.0);
  495.     if (arg < MINEXP)
  496.     return -1.0;
  497.     arg = gp_exp(arg);
  498.  
  499.     /* Choose infinite series or continued fraction. */
  500.  
  501.     if ((x > 1.0) && (x >= a + 2.0)) {
  502.     /* Use a continued fraction expansion */
  503.  
  504.     double          rn;
  505.     double          rnold;
  506.  
  507.     aa = 1.0 - a;
  508.     b = aa + x + 1.0;
  509.     pn1 = 1.0;
  510.     pn2 = x;
  511.     pn3 = x + 1.0;
  512.     pn4 = x * b;
  513.     rnold = pn3 / pn4;
  514.  
  515.     for (i = 1; i <= ITMAX; i++) {
  516.  
  517.         aa++;
  518.         b += 2.0;
  519.         an = aa * (double) i;
  520.  
  521.         pn5 = b * pn3 - an * pn1;
  522.         pn6 = b * pn4 - an * pn2;
  523.  
  524.         if (pn6 != 0.0) {
  525.  
  526.         rn = pn5 / pn6;
  527.         if (fabs(rnold - rn) <= min(MACHEPS, MACHEPS * rn))
  528.             return 1.0 - arg * rn * a;
  529.  
  530.         rnold = rn;
  531.         }
  532.         pn1 = pn3;
  533.         pn2 = pn4;
  534.         pn3 = pn5;
  535.         pn4 = pn6;
  536.  
  537.         /* Re-scale terms in continued fraction if terms are large */
  538.         if (fabs(pn5) >= OFLOW) {
  539.  
  540.         pn1 /= OFLOW;
  541.         pn2 /= OFLOW;
  542.         pn3 /= OFLOW;
  543.         pn4 /= OFLOW;
  544.         }
  545.     }
  546.     } else {
  547.     /* Use Pearson's series expansion. */
  548.  
  549.     for (i = 0, aa = a, an = b = 1.0; i <= ITMAX; i++) {
  550.  
  551.         aa++;
  552.         an *= x / aa;
  553.         b += an;
  554.         if (an < b * MACHEPS)
  555.         return arg * b;
  556.     }
  557.     }
  558.     return -1.0;
  559. }
  560.  
  561. /***********************************************************************
  562.      double ranf(double init)
  563.                 RANDom number generator as a Function
  564.      Returns a random floating point number from a uniform distribution
  565.      over 0 - 1 (endpoints of this interval are not returned) using a
  566.      large integer generator.
  567.      This is a transcription from Pascal to Fortran of routine
  568.      Uniform_01 from the paper
  569.      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
  570.      with Splitting Facilities." ACM Transactions on Mathematical
  571.      Software, 17:98-111 (1991)
  572.  
  573.                GeNerate LarGe Integer
  574.      Returns a random integer following a uniform distribution over
  575.      (1, 2147483562) using the generator.
  576.      This is a transcription from Pascal to Fortran of routine
  577.      Random from the paper
  578.      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
  579.      with Splitting Facilities." ACM Transactions on Mathematical
  580.      Software, 17:98-111 (1991)
  581. ***********************************************************************/
  582. static double          ranf(init)
  583. double init;
  584. {
  585.  
  586.     long            k, z;
  587.     static int      firsttime = 1;
  588.     static long     s1, s2;
  589.  
  590.     /* (Re)-Initialize seeds if necessary */
  591.     if (init < 0.0 || firsttime == 1) {
  592.     firsttime = 0;
  593.     s1 = 1234567890L;
  594.     s2 = 1234567890L;
  595.     }
  596.     /* Generate pseudo random integers */
  597.     k = s1 / 53668L;
  598.     s1 = Xa1 * (s1 - k * 53668L) - k * 12211;
  599.     if (s1 < 0)
  600.     s1 += Xm1;
  601.     k = s2 / 52774L;
  602.     s2 = Xa2 * (s2 - k * 52774L) - k * 3791;
  603.     if (s2 < 0)
  604.     s2 += Xm2;
  605.     z = s1 - s2;
  606.     if (z < 1)
  607.     z += (Xm1 - 1);
  608.  
  609.     /*
  610.      * 4.656613057E-10 is 1/Xm1.  Xm1 is set at the top of this file and is
  611.      * currently 2147483563. If Xm1 changes, change this also.
  612.      */
  613.     return (double) 4.656613057E-10 *z;
  614. }
  615.  
  616. /* ----------------------------------------------------------------
  617.    Following to specfun.c made by John Grosh (jgrosh@arl.mil)
  618.    on 28 OCT 1992.
  619.    ---------------------------------------------------------------- */
  620.  
  621. void f_normal()    /* Normal or Gaussian Probability Function */
  622. {
  623. struct value a;
  624. double x;
  625.  
  626.     /* ref. Abramowitz and Stegun 1964, "Handbook of Mathematical 
  627.        Functions", Applied Mathematics Series, vol 55,
  628.        Chapter 26, page 934, Eqn. 26.2.29 and Jos van der Woude 
  629.            code found above */
  630.  
  631.     x = real(pop(&a));
  632.  
  633. #ifndef ERF
  634.         x = 0.5 * SQRT_TWO * x;
  635.         x = 0.5 * (1.0 + (x < 0.0 ? -igamma(0.5, x * x) : igamma(0.5, x * x)));
  636. #else
  637.     x = 0.5 * (1.0 + erf(0.5 * SQRT_TWO * x));
  638. #endif
  639.         push( Gcomplex(&a,x,0.0) );
  640. }
  641.  
  642. void f_inverse_normal()  /* Inverse normal distribution function */
  643. {
  644. struct value a;
  645. double x;
  646. double inverse_normal_func();
  647.  
  648.     x = real(pop(&a));
  649.  
  650.     if (fabs(x) >= 1.0) {
  651.         undefined = TRUE;
  652.         push(Gcomplex(&a,0.0, 0.0));
  653.     } else {
  654.         push( Gcomplex(&a,inverse_normal_func(x), 0.0) );
  655.     }
  656. }
  657.  
  658.  
  659. void f_inverse_erf()  /* Inverse error function */
  660. {
  661. struct value a;
  662. double x;
  663. double inverse_error_func();
  664.  
  665.     x = real(pop(&a));
  666.  
  667.     if (fabs(x) >= 1.0) {
  668.         undefined = TRUE;
  669.         push(Gcomplex(&a,0.0, 0.0));
  670.     } else {
  671.         push( Gcomplex(&a,inverse_error_func(x), 0.0) );
  672.     }
  673. }
  674.  
  675. static double 
  676. inverse_normal_func(p)
  677. double p;
  678. {
  679.     /* 
  680.            Source: This routine was derived (using f2c) from the 
  681.                    FORTRAN subroutine MDNRIS found in 
  682.                    ACM Algorithm 602 obtained from netlib.
  683.  
  684.                    MDNRIS code contains the 1978 Copyright 
  685.                    by IMSL, INC. .  Since MDNRIS has been 
  686.                    submitted to netlib it may be used with 
  687.                    the restriction that it may only be 
  688.                    used for noncommercial purposes and that
  689.                    IMSL be acknowledged as the copyright-holder
  690.                    of the code.
  691.         */
  692.  
  693.     /* Initialized data */
  694.     static double eps = 1e-10;
  695.     static double g0 = 1.851159e-4;
  696.     static double g1 = -.002028152;
  697.     static double g2 = -.1498384;
  698.     static double g3 = .01078639;
  699.     static double h0 = .09952975;
  700.     static double h1 = .5211733;
  701.     static double h2 = -.06888301;
  702.     static double sqrt2 = 1.414213562373095;
  703.  
  704.     /* Local variables */
  705.     static double a, w, x;
  706.     static double sd, wi, sn, y;
  707.  
  708.     double inverse_error_func();
  709.  
  710.     /* Note: 0.0 < p < 1.0 */
  711.  
  712.     /* p too small, compute y directly */
  713.     if (p <= eps) {
  714.         a = p + p;
  715.         w = sqrt(-(double)log(a + (a - a * a)));
  716.  
  717.         /* use a rational function in 1.0 / w */
  718.         wi = 1.0 / w;
  719.         sn = ((g3 * wi + g2) * wi + g1) * wi;
  720.         sd = ((wi + h2) * wi + h1) * wi + h0;
  721.         y = w + w * (g0 + sn / sd);
  722.         y = -y * sqrt2;
  723.     } else {
  724.         x = 1.0 - (p + p);
  725.         y = inverse_error_func(x);
  726.         y = -sqrt2 * y;
  727.     }
  728.     return(y);
  729. } 
  730.  
  731.  
  732. static double 
  733. inverse_error_func(p) 
  734. double p;
  735. {
  736.     /* 
  737.            Source: This routine was derived (using f2c) from the 
  738.                    FORTRAN subroutine MERFI found in 
  739.                    ACM Algorithm 602 obtained from netlib.
  740.  
  741.                    MDNRIS code contains the 1978 Copyright 
  742.                    by IMSL, INC. .  Since MERFI has been 
  743.                    submitted to netlib, it may be used with 
  744.                    the restriction that it may only be 
  745.                    used for noncommercial purposes and that
  746.                    IMSL be acknowledged as the copyright-holder
  747.                    of the code.
  748.         */
  749.  
  750.  
  751.  
  752.     /* Initialized data */
  753.     static double a1 = -.5751703;
  754.     static double a2 = -1.896513;
  755.     static double a3 = -.05496261;
  756.     static double b0 = -.113773;
  757.     static double b1 = -3.293474;
  758.     static double b2 = -2.374996;
  759.     static double b3 = -1.187515;
  760.     static double c0 = -.1146666;
  761.     static double c1 = -.1314774;
  762.     static double c2 = -.2368201;
  763.     static double c3 = .05073975;
  764.     static double d0 = -44.27977;
  765.     static double d1 = 21.98546;
  766.     static double d2 = -7.586103;
  767.     static double e0 = -.05668422;
  768.     static double e1 = .3937021;
  769.     static double e2 = -.3166501;
  770.     static double e3 = .06208963;
  771.     static double f0 = -6.266786;
  772.     static double f1 = 4.666263;
  773.     static double f2 = -2.962883;
  774.     static double g0 = 1.851159e-4;
  775.     static double g1 = -.002028152;
  776.     static double g2 = -.1498384;
  777.     static double g3 = .01078639;
  778.     static double h0 = .09952975;
  779.     static double h1 = .5211733;
  780.     static double h2 = -.06888301;
  781.  
  782.     /* Local variables */
  783.     static double a, b, f, w, x, y, z, sigma, z2, sd, wi, sn;
  784.  
  785.     x = p;
  786.  
  787.     /* determine sign of x */
  788.     if (x > 0)
  789.         sigma = 1.0;
  790.     else
  791.         sigma = -1.0;
  792.  
  793.     /* Note: -1.0 < x < 1.0 */
  794.  
  795.     z = fabs(x);
  796.  
  797.     /* z between 0.0 and 0.85, approx. f by a 
  798.        rational function in z  */
  799.  
  800.     if (z <= 0.85) {
  801.         z2 = z * z;
  802.         f = z + z * (b0 + a1 * z2 / (b1 + z2 + a2 
  803.             / (b2 + z2 + a3 / (b3 + z2))));
  804.  
  805.     /* z greater than 0.85 */
  806.     } else {
  807.         a = 1.0 - z;
  808.         b = z;
  809.  
  810.         /* reduced argument is in (0.85,1.0), 
  811.            obtain the transformed variable */
  812.  
  813.         w = sqrt(-(double)log(a + a * b));
  814.  
  815.         /* w greater than 4.0, approx. f by a 
  816.            rational function in 1.0 / w */
  817.  
  818.         if (w >= 4.0) {
  819.             wi = 1.0 / w;
  820.             sn = ((g3 * wi + g2) * wi + g1) * wi;
  821.             sd = ((wi + h2) * wi + h1) * wi + h0;
  822.             f = w + w * (g0 + sn / sd);
  823.  
  824.         /* w between 2.5 and 4.0, approx. 
  825.            f by a rational function in w */
  826.  
  827.         } else if (w < 4.0 && w > 2.5) {
  828.             sn = ((e3 * w + e2) * w + e1) * w;
  829.             sd = ((w + f2) * w + f1) * w + f0;
  830.             f = w + w * (e0 + sn / sd);
  831.  
  832.         /* w between 1.13222 and 2.5, approx. f by 
  833.            a rational function in w */
  834.         } else if (w <= 2.5 && w > 1.13222) {
  835.             sn = ((c3 * w + c2) * w + c1) * w;
  836.             sd = ((w + d2) * w + d1) * w + d0;
  837.             f = w + w * (c0 + sn / sd);
  838.         }
  839.     }
  840.     y = sigma * f;
  841.     return(y);
  842. }
  843.  
  844.